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Abstract 



The dynamics of first order phase transitions are studied in the context 
of (3+l)-dimensional scalar field theories. Particular attention is paid to 
the question of quantifying the strength of the transition, and how 'weak' 
and 'strong' transitions have different dynamics. We propose a model with 
two available low temperature phases separated by an energy barrier so 
that one of them becomes metastable below the critical temperature Tc- 
The system is initially prepared in this phase and is coupled to a ther- 
mal bath. Investigating the system at its critical temperature, we find 
that 'strong' transitions are characterized by the system remaining local- 
ized within its initial phase, while 'weak' transitions are characterized by 
considerable phase mixing. Always at we argue that the two regimes 
are themselves separated by a (second order) phase transition, with an or- 
der parameter given by the fractional population difference between the 
two phases and a control parameter given by the strength of the scalar 
field's quartic self-coupling constant. We obtain a Ginzburg-like criterion 
to distinguish between 'weak' and 'strong' transitions, in agreement with 
previous results in (2-|-l)-dimensions. 
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I. INTRODUCTION 



The fact that the gauge symmetries describing particle interactions can be restored 
at high enough temperatures has led, during the past 15 years or so, to an active research 
program on the possible implications that this symmetry restoration might have had to 
the physics of the very early Universe. One of the most interesting and popular possi- 
bilities is that during its expansion the Universe underwent a series of phase transitions, 
as some higher symmetry group was successively broken into products of smaller groups, 
up to the present standard model described by the product SU{3)c <S> SU{2)l ® U{1)y- 
Most models of inflation and the formation of topological (and nontopological defects) 
are well-known consequences of taking the existence of cosmological phase transitions 
seriously |jl|]. 

One of the motivations of the present work comes from the possibility that the baryon 
asymmetry of the Universe could have been dynamically generated during a first order 
electroweak phase transition ||2[. As is by now clear, a realistic calculation of the net 
baryon number produced during the transition is a formidable challenge. We probably 
must invoke physics beyond the standard model (an exciting prospect for most people) 
0, push perturbation theory to its limits (and beyond, due to the nonperturbative nature 
of magnetic plasma masses that regulate the perturbative expansion in the symmetric 
phase), and we must deal with nonequilibrium aspects of the phase transition. Here we 
will focus on the latter problem, as it seems to us to be the least discussed of the pillars 
on which most baryon number calculations are built upon. To be more specific, we can 
separate the nonequilibrium aspects of the phase transition into two main subdivisions. 
If the transition proceeds by bubble nucleation, we can study the propagation of bubbles 
in the hot plasma and the transport properties of particles through the bubble wall. A 
considerable amount of work has been devoted to this issue, and the reader can consult 
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the works of Ref. 0] for details. These works assume that homogeneous nucleation theory 
is adequate to investigate the evolution of the phase transition, at least for the range of 
parameters of interest in the particular model being used to generate the baryon asym- 
metry. This brings us to the second important aspect of the nonequilibrium dynamics 
of first order phase transitions, namely the validity of homogeneous nucleation theory to 
describe the approach to equilibrium. This is the issue addressed in the present work. 

Nucleation theory is a well-studied, but far from exhausted, subject. Since the pio- 
neering work of Becker and Doring on the nucleation of droplets in supercooled vapor [^, 
the study of first order phase transitions has been of interest to investigators in several 
fields, from meteorology and materials science to quantum field theory and cosmology. 
Phenomenological field theories were developed by Cahn and Hilliard and by Langer 
in the context of coarse-grained time-dependent Ginzburg-Landau models, in which an 
expression for the decay rate per unit volume was obtained by assuming a steady-state 
probability current flowing through the saddle-point of the free-energy functional 
The application of metastable decay to quantum field theory was initiated by Voloshin, 
Kobzarev, and Okun and soon after put onto firmer theoretical ground by Coleman 
and Callan The generalization of these results for finite temperature field theory 



was first studied by Linde [jrT| , and has been the focus of much recent attention [T^ . 



The crucial ingredient in the evaluation of the decay rate is the computation of the 
imaginary part of the free energy. As shown by Langer [^, the decay rate TZ is propor- 
tional to the imaginary part of the free energy JF, 

7^ = ^^Im-F , (1) 

where is the negative eigenvalue related to metastability, which depends on nonequi- 
librium aspects of the dynamics, such as the coupling strength to the thermal bath. Since 
= —TlnZ, where Z is the partition function, the computation for the rate boils down 



to the evaluation of the partition function for the system comprised of bubbles of the 
lower energy phase inside the metastable phase. For a dilute gas of bubbles only, the 
partition function for several bubbles is given by [0,0, 
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where is the metastable vacuum field configuration and (^b is the bubble configu- 
ration, the bounce solution to the 0(3)-symmetric Euclidean equation of motion. We 
must evaluate the partition functions above. This is done by the saddle-point method, 
expanding the scalar field 0(x, r), such that 0(x, r) iff + ({'^,t) for Z{ipf), and 
0(x, r) — i> v^6(x) + ?7(x, r) for Z{ipb), where Ci'^j'^) vi^:'^) small fluctuations 
about equilibrium. Skipping details [jl2|, up to one-loop order one obtains for the ratio 
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Z[(pb) 1— loop order 

Z{^f) 



det(-nB + V"{^b))p 



-AS 



(3) 



_det{-nE + V"iipf))p_ 
where [det(M)^]-5 = J Dr] exp [- J/^ dr J d^xl7][M]7]} and A 5 = Ssi^b) - Sei^f) is the 
difference between the Euclidean actions for the field configurations (pb and iff. [Note 
that Ssif), and hence AS", does not include any temperature corrections.] Thus, the 
free energy of the system is, 

^det{-nE + V"iipb))i 
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We are briefly reproducing this computation here because we want to stress the im- 
portance of the assumptions built into it. First, that the partition function is given by 
equation ^ within the dilute gas approximation, and second, that the partition function 
is evaluated assuming small fluctuations about the homogeneous metastable state iff. 
It is clear that for situations in which there are large amplitude fluctuations about the 



metastable equilibrium state the above formula must break down. Thus the breakdown 
of the expression for the rate is intimately connected with the question of how well- 
localized the system is about the metastable state as the temperature drops below the 
critical temperature Tc. 

This question has been addressed in the context of the electroweak phase transi- 
tion in the works listed in Ref. fl^. The common assumption of these works is that 



for weak enough transitions the dynamics is dominated by correlation-volume large- 
amplitude thermal fluctuations, dubbed subcritical bubbles, which promote considerable 
phase mixing as the Universe cools below the critical temperature. Within the validity of 
this analytical approach it was shown that homogeneous nucleation is only justified for 
Higgs masses below 70 GeV or so, which is both dangerously close to the present lower 



bound on the Higgs mass |T5[ , and to the limit of validity of the perturbatively evaluated 



effective potential [0. Furthermore, as with any analytical treatment of nonequilib- 
rium dynamics, many aspects of the complicated kinetics of the system are not included. 
For example, even though correlation-volume bubbles may be the dominant fluctuations 
about equilibrium, there will be bubbles of different sizes present which may percolate 
and also acquire some thermal velocity due to diffusive processes. It is clear that a final 
answer to the problem can only be given by a combination of analytic and numerical 
tools. 

In order to clarify the distinction between 'weak' and 'strong' first order transitions, 
one of us has recently studied the nonequilibrium dynamics of a (2+l)-dimensional model 
of a scalar field coupled to a thermal bath The nonlinear interactions were chosen 
to reflect the gross properties of the electroweak effective potential, although the model 
only deals with a real scalar field. It was shown that there is a very clear distinction 
between a weak and a strong transition, and that one should expect a very different 
dynamics between the two cases. In the present work, we generalize these results to (3+1)- 



dimensions and also address important issues concerning the reliability of the numerical 
results. As we will show, it is clear that the two regimes are easily distinguishable, as in 
the (2+l)-dimensional case. Using homogeneous nucleation can easily lead to the wrong 
description of the dynamics. 

In passing, we note that somewhat similar results have been obtained in the context of 
binary-fluid mixtures, where the (conserved) order parameter is the local concentration 
of one of the components of the mixture . If the system is quenched to concentrations 
above the spinodal (the inflection point in the free-energy density), the transition evolves 
by spinodal decomposition; otherwise, nucleation occurs. The transition between the two 
regimes was shown in that case to be smooth. However, we must remember that here we 
are dealing with a non-conserved order parameter and have much faster dynamics than 
in binary fluid mixtures. As we will show, the transition between the two regimes is more 
dramatic in our context. 

This paper is organized as follows. In the next Section we introduce the model we 
will use in the numerical simulations. In Section 3 we discuss details of the numerical 
approach used to study the dynamics of the system, including the implementation of 
the code on a parallel machine. In Section 4 we discuss the numerical results and their 
reliability. We conclude in Section 5 with a summary of our results and an outlook into 
future work. 



II. THE MODEL 

The homogeneous part of the free energy density is written as 

U{<P, T) = I (T^ - T|) 0^ - |t0^ + . (5) 

This choice intentionally resembles the electroweak effective potential to some order in 
perturbation theory, although here (/)(x, t) is a real scalar field, as opposed to the magni- 
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tude of the Higgs field. The goal is to explore the possible dynamics of a model described 
by the above free-energy density, generalizing the results obtained in Ref. to (3+1)- 



dimensions. The analogy with the electroweak model is suggestive but not quantitative. 

Introducing dimensionless variables x = a^T2X, i = a^T2t, X = a~^T2^(f), and 
= T/T2, the Hamiltonian is, 



1 1 yx |2 +1 (e'^ -i)x'- ^ex' + -x^ 



2 ' ' ' 2 V / 3 4 



(6) 



3 ~ 1 

where a = a~4a, and A = a~2X (henceforth we drop the tildes). For temperatures above 
61 = (1 — «^/4A)^2 there is only one minimum at X = 0. At ^ = ^1 an infiection point 
appears at X^^f = a6i/2X. Below 61 the inflection point separates into a maximum and 



a minimum given by X^ = tf liyl — 4A(1 — 1/0^) jo? . At the critical temperature 
= (1 - 2aV9A)-2 the two minima, at X^ = and X+ are degenerate. Below 6^ the 
minimum at X+ becomes the global minimum and the Xo-phase becomes metastable. 
Finally, at ^ = 1 the barrier between the two phases at X_ disappears. 

The coupling with the thermal bath will be modelled by a Markovian Langevin equa- 
tion which, in terms of the dimensionless variables defined above, is 

d'^X „2.. dX dU(X,e) , 

where t] is the dimensionless viscosity coefficient, and ^ the dimensionless stochastic noise 
with vanishing mean, related to 77 by the fiuctuation-dissipation theorem, 

(e(x, t)C(x', t')) = 2^e5{t - t')5\^ - x') . (8) 

A few comments are in order concerning our choice of equation. It is clear that we are 
assuming that X(x, t) represents the long-wavelength modes of the scalar field. When- 
ever one discretizes a continuum system there is an implicit coarse-graining scale built 
in. We encapsulate information about the shorter-wavelength modes, which have faster 



relaxation time-scales, in the dissipation and noise terms. In principle it should be pos- 
sible to derive an effective Langevin-like equation for the slow modes by integrating out 
the fast modes from the effective action. This is a complicated problem, and progress has 
been slow. Recent work indicates that one should expect departures from the Langevin 
equation written above [^, although details are sensitive to the particular model one 



starts with. For example, the noise may be colored (with more complicated correlation 
functions) and the coupling to the bath may be multiplicative, as opposed to the additive 
coupling chosen above. Here, we will adopt the above equation as a first step. We do 
not expect that the nature of the noise will change the final equilibrium properties of the 
system, but mostly the relevant relaxation time-scales. Since the physical results here 
are related to the final equilibrium state of the system, we believe that they will not be 
affected by more complicated representations of the coupling of the field to the thermal 
bath. However, a more thorough examination of this question deserves further study. 

A related topic is the choice of coarse-graining scale, which is embedded in the lattice 
spacing used in the simulations. It is well-known that any classical field theory in more 
than one spatial dimension is ultra-violet divergent, and that the lattice spacing serves 
as an ultra-violet cutoff. This being the case, one should be careful when mapping from 
the lattice to the continuum theory. If one is to probe physics at shorter wavelengths, 
renormalization counterterms should be included in the lattice formulation so that a 
proper continuum limit is obtained on the lattice within the validity of perturbation 
theory. This point has been emphasized in Ref. |2^, where a (2+l)-dimensional study 



of nucleation was performed for a temperature-independent potential. Renormalization 
counterterms (of order 9 In 5a; for lattice spacing 5x) for a particular renormalization 
prescription were obtained, and the results shown to be lattice-space independent. 

Here, due to the temperature dependence of the potential, the renormalization pre- 
scription of Ref. |2y] does not work. Instead, we will use 5x = 1 throughout this work. 
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It turns out that for all cases studied the mean- field correlation length = V"{Xo, 9c) 
is sufficiently larger than unity to justify this choice. Modes with shorter wavelengths 
are coupled through the noise into the dynamics of the longer wavelength modes, as 
described by equation |^ above. 



III. LATTICE FORMULATION 



The system is now discretized onto a lattice of length L with grid spacing 6x, time step 
St, and total run time At. Using a standard second order staggered leapfrog approach 
equation ^ becomes 

1 



X 



i,n+l/2 



1 + ^¥r. 
,n+l — Xi^n + St 



1 - -r]tn] Xi^n-l/2 + St V'Xi,„ - — 



(9) 



where i-indices are spatial and n-indices temporal. The continuum white noise ^(x, t) is 
replaced by its discretized analogue ^i,„ by requiring that the discrete noise be uncorre- 
lated on all scales above the shortest simulated. The discretised form of the ffuctuation- 
dissipation relation of equation || becomes 



St [Sx] 
The discrete noise is hence approximated by 



(10) 



2r]e 



(11) 



^ St {Sxf 

where Qi^n is a unit- variance Gaussian random number at each point in lattice spacetime. 

Since we are modelling an unbounded system we do not want our simulation volume 
to have a distinct surface; we therefore use periodic boundary conditions. However, such 
boundary conditions may induce errors if a simulation runs for longer than a time causally 
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equivalent to L/2 as the periodicity then introduces spurious long-range correlations. 
Since we must run very long simulations to guarantee equilibration, were this constraint 
to apply we would be forced to use impossibly large lattices. Fortunately, the presence of 
the noise term, uncorrelated at each point in lattice spacetime, has the effect of swamping 
any such effect. 

We must now run our simulations on large lattices (to reduce finite size effects) many 
times over (to reduce statistical noise) and for long run times (to ensure reaching equilib- 
rium). Typically any attempt to reduce one of these constraints is met by a corresponding 
increase in another — for example, smaller grids give noisier results requiring more runs. 
As a consequence, we soon find ourselves at the limits of what is possible on a worksta- 
tion. We have therefore parallehsed the code and implemented it on a 128-node APIOOO. 
The overall lattice is subdivided into an appropriate number of sub-lattices, each of which 
is local to a single node. These sub-lattices are defined with an overlap, such that each 
edge of any node's sub-lattice is included within the body of one of its neighbours. At 
each time step each node evolves the body of its local sub-lattice, but not the edges, for 
which insufficient data is locally available to calculate spatial derivatives. Each node then 
swaps the necessary data to update the overlapping edges of their associated sub-lattices 
with each of its neighbouring nodes. 

The one qualitatively different feature of the parallel code is in the implementation of 
the random number generator. Computational random number generators are required 
to produce the same sequence of numbers whenever they are given the same initial 
conditions. Therefore what they actually generate, at least ideally, is a predictable, 
periodic, sequence of pseudo-random numbers, any sufficiently short {i.e. substantially 
shorter than the period) subset of which has statistical properties indistinguishable from 
a genuine random sequence. In our case we need random numbers at every point in 
lattice spacetime which are uncorrelated across the entire simulation spacetime; each node 
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requires numbers which are random not only locally at the node itself, but also across all 
the other nodes too. Either each node must have a different generator (highly impractical 
for any more than a few nodes) or each node must be allocated a unique subsequence of a 
single generator's full sequence. Most generators are based on an iterative scheme, where 
each number in the sequence is calculated from some of the previous numbers in the 
sequence. However, if these previous numbers are not members of the local subsequence 
then there is a communication cost incurred in fetching them from the relevant node. 
In order to maximise the efficiency of the code, we require a generator whose sequence 
can be broken down into subsequences with elements generated by reference to previous 
members of that subsequence alone. Moreover, for large lattices and long run times we 
require a generator with an extremely long period. Thus for L = 48, 6t = 0.1, and 
a running to At = 3000 (values used below) we require of the order of 2^^ random 
numbers, and hence a generator with a period many times longer than this. A solution 
to this problem, applicable across any 2" nodes, with a period of 2^*^^ — 1 and with 
each element in each subsequence calculable completely locally, is given by a particularly 
elegant parallelisation of the generalised feedback shift register algorithm |^T|; this is the 
generator implemented here. 



IV. NUMERICAL EXPERIMENT AND RESULTS 

As pointed out in the Introduction, the question as to whether a first order phase 
transition is 'weak' or 'strong' boils down to how well localized in the metastable state 
the system is as the temperature drops below the critical temperature. In order to 



address this question, following the procedure of Ref. |T^, we will study the behavior 
of the system at the critical temperature, when the two minima are degenerate. The 
reason for this choice follows naturally from the fact that we are interested on the way 
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by which the system approaches equihbrium as the temperature drops below T^. The 
detailed dynamics will depend on the relative fraction of the total volume occupied by 
each phase; if at the system is well localized about the X = minimum, as the 
temperature drops the transition may evolve by nucleation and subsequent percolation 
of bubbles larger than a critical size. If, on the other hand, considerable phase-mixing 
occurs already at T^, we expect the transition to evolve by domain coarsening, with the 
domains of the X+ phase eventually permeating the whole volume. 

Let us call the two phases the 0-phase and the +-phase, corresponding to the local 
equilibrium values X = Xo = 0, and X = X^, respectively. We can quantify the phase 
distribution of the system as it evolves according to equation |^, by measuring the fraction 
of the total volume in each phase. This is done by simply counting the total volume of 
the system at the left of the potential barrier's maximum height {i.e., X < X_ = -^max), 
corresponding to the 0-phase. Dividing by the total volume, we obtain the fraction of 
the system in the 0-phase, fo{t), such that 

m+uit) = i , (12) 

where, of course, corresponds to the fractional volume in the +-phase. 

A further measure of any configuration is given by the volume-averaged order param- 
eter, {X)(t) = J dV X{t). A localised configuration (/^'^ > 0.5) then corresponds 
to {Xy'^ < Xjnax, and a fully phase-mixed configuration (/^'^ ~ 0.5) to {Xy'^ = Xmax, 
where the superscript 'eq' refers to final ensemble-averaged equilibrium values of fo{t) 
and (X)(t). 

We prepare the system so that initially it is well localised in the 0-phase, with /o(0) = 
1 and (X)(0) = 0. These initial conditions are clearly the most natural choice for the 
problem at hand. If one has cosmology in mind, it is quite possible that as the system 
slowly cools down (we are not interested in phase transitions close to the Planck scale), 
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fluctuations from the high temperature phase X = to the phase are aheady occuring 
before T^. is reached. (In this case, our arguments are even stronger.) However, we will 
adopt the best-case scenario for homogeneous nucleation to work, in which the system 
managed to reach the X = phase homogeneously, so that the initial state is a thermal 
state with mean at Xq. If one has more concrete applications in mind, we can assume 
that we quenched the system to its critical temperature, making sure that the order 
parameter remains localized about the high-temperature phase. Since thermalization 
happens very fast in the simulations, the exact point by point initial conditions should 
not be important, and we can view the flrst few time steps as generating an initial 
thermal distribution with /o(0) ~ 1 and (X) ~ 0, so that the average kinetic energy 
per lattice point satisfles the equipartition theorem, jjE^ = |T. For simplicity we take 
X = 0, X = everywhere initially. 

There are two parameters controlling the strength of the transition, a and A. In the 
previous (2-|-l)-dimensional work, a was chosen to vary while A was kept flxed. It is 
really immaterial which parameter is held flxed, or if both are made to vary, but in order 
to keep closer to the spirit of the electroweak model we will flx a and let A vary. As is 
well-known, A is related to the Higgs mass, while a is related to the gauge-boson masses 
101 . The connection with the electroweak model is straightforward. If we consider as an 
example the unimproved one-loop approximation, the effective potential is |@], 



where D and E are given by D = [6{Mw/af + ^{Mz/af + Q^Mx/af] /24 ~ 0.17 and 
E = [6{Mw/(tY + ^iMz/af] jVl-n ~ 0.097, for Mw = 80.6 GeV, Mz = 91.2 GeV, 
Mr = 174 GeV and a = 246 GeV. T2 is given by. 



where the physical Higgs mass is given in terms of the 1-loop corrected A as Mfj = 




(13) 




(14) 



13 



(2A + 125) (t2, with B = (6M^ + 3Af| - 12M|) /647^V^ and the temperature-corrected 
Higgs self-couphng is, 
1 



At — A 



j:g, (^)\n{Ml/csT^)-j:9F (^)' In (m|/c,T^ 

B F 



(15) 



167r2 

where the sum is performed over bosons and fermions (in our case only the top quark) 
with their respective degrees of freedom gB{F), and Inc^ = 5.41 and Inci? = 2.64. 
Thus, the correspondence with our (dimensionless) parameters is 

3E A 

a = J = 0.065 , and A = — ^ = 1.72Ar . (16) 

(2D)t {2D)-2 

Once this is estabhshed, the numerical experiment proceeds as follows: i) Choose a = 

0.065; ii) Prepare the system in the initial state described above, and measure the value 

of fo(t) and {X){t) for several values of A, as the system evolves according to equation |^ 

As with any numerical experiment, one must make sure that the results are independent 

of lattice artifacts (or at least the dependence is understood), such as its finite size and 

choice of time-step for the evolution routine. Before we discuss our results, we will address 

these issues in the following Subsection. 



A. Finite Lattices and the Thermodynamic Limit 

Whenever simulating a system on the lattice, we are faced with the hard decision of 
having to achieve a compromise between approximating the infinite volume limit, and 
having fairly reasonable integration times. This problem is particularly serious in the 
context of phase transitions, as it is well-known that symmetry breaking only occurs in 
the infinite volume limit; at finite volumes, there is a nonzero probability that a large 
fluctuation will restore the broken symmetry. Even though this is formally true, we will 
argue here that this does not represent a problem to our simulations, if we make sure 
that the lattice is large enough. There is a large amount of literature on finite-size effects 
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and how they are handled in different contexts |2^, and we do not intend to reproduce 
these results. What we want to do is to bring this issue closer to our problem. 

We are interested in studying the system given by the free-energy density of equation 
^, at the critical temperature 9c when the two minima are degenerate. The system is 
prepared in the 0-phase, and we measure the fraction of the volume in each phase as 
it evolves. We will give a rough estimate of how large the lattice should be in order to 
suppress spurious symmetry restoration (that is, /o — > 0.5, in our case) due to the lattice 
size. There are two relevant time-scales in the problem, the relaxation time-scale for 
small-amplitude fluctuations within the X = well, Trei, and the 'escape' time-scales for 
large amplitude fluctuations into the +-phase, t^sc- In terms of the rate per unit volume 
for each process, we write the relevant time-scales as 

yr,ei ~ r,;/ ~ T7-? , (17) 

where 7rei is the typical relaxation time-scale for short amplitude fluctuations in units of 
T~^, and 

VV^^^r^r-^lr^TeM-Ff/T] , (18) 

where Fj is the free energy of the large amplitude fluctuation. The condition for large 
amplitude fluctuations to be suppressed in comparison to typical relaxation processes is 
then, 

— >1 ^ exp[F^/T] >7rei . (19) 

''"rel 

To estimate Tesc note that within the Gaussian approximation, a homogeneous fluctuation 
of volume V and amplitude 0a about equilibrium (0 = 0, for simplicity) has free energy 

F;(0A,\/,T) = ^m2(T)0^ , (20) 
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where im?{T) = V"{(j) = 0,T) (we neglect the gradient contribution, as it would suppress 
the fluctuation even further, making our arguments stronger). We are interested in fluc- 
tuations about X = (we now go back to our dimensionless variables), at the critical 
temperature 9c = {I — 2a^/9A)~5. We expect growing instabilities to be triggered when- 
ever fluctuations probe the nonlinearities in the free energy. This is corroborated by the 
results in Ref. [|l7l and, as we will soon see, here also. Thus we consider fluctuations with 
amplitude equal to the nearest inflection point to X = 0, namely Xa = ^ (l — 
Writing their volume V = ^ (^0^ terms of the correlation length ^{9c) = {9l — 1)^^, 
with n a real number, we obtain 



F 



a 



^^ 0.088^n^ . (21) 

_ 1 

From our arguments above it is clear that 7rei = (6*^ — 1) ^ . Let us consider an example 
which will be relevant later on. Take a = 0.065 and A = 0.020. In this case, the 
correlation length is ^(^c) — 4.5 and we obtain, from equation |19|, n ^ 0.91. Since the 
radius of the fluctuation is Rf = this result implies that fluctuations probing the 
inflection point with radius Rf = 0.91 x 4.5 ^ 4.1 are probable within typical relaxation 
time scales of the system. Thus the lattice length L should be sufficiently larger than 
about 2Rf {L > 10 or so) to avoid spurious symmetry restoration. For L = 20 the volume 
ratio of the lattice to the above fluctuations is around L^/AR^ ~ 25, and such processes 
alone cannot restore the symmetry within time-scales of interest in the dynamics. In all 
the results quoted in this work we use L = 48. 

We give two pieces of evidence supporting these arguments. In Fig.l the equilibrium 
values of the 0-phase fraction, f^"^, and of the volume-averaged order parameter, {Xy^, 
are given as a function of the lattice size L, for a = 0.065 and A = 0.020. A is chosen so 
that for large lattices the system remains localized mostly in the 0-phase. Note that as L 
decreases f^'^ 0.5, representing a spurious 'symmetry restoration' due to the smallness 
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of the lattice. For small lattices, the two phases are completely mixed. In fact, changes 
can be seen between L = 10 and L = 20, in qualitative agreement with our arguments 
above. Note that for large enough lattices, the equilibrium value approaches a stable 
value which is independent of the lattice size. For all practical purposes, this is the 
infinite volume limit. Fluctuations large enough to restore the symmetry are possible, 
but with negligible probability. 

The reader may wonder why for large enough lattices f^'^ ^ 1. The reason for this is 
that at finite temperatures there is a nonvanishing probability per unit volume of having 
thermal fluctuations populating the +-phase. Even though these fluctuations are unstable 
and shrink away, there will be an equilibrium distribution of bubbles suppressed by a 
Boltzmann factor. For very strong transitions (very small A), f^^ ~ 1, and a negligible 
fraction of the system is in the +-phase. We refer the reader to the paper by Gelmini and 
Gleiser, Ref. for details. In Fig. 2 we show a phase space portrait of the system for a 



given point on the lattice for different lattice sizes. That is, we choose a particular point 
X{xo, Ho, Zo, t) and follow its evolution, making a plot of X(xo, i/o, ^o, t) vs X{xo, i/o, Zo, t). 
Taking a = 0.065 and A = 0.010, Fig. 2a shows the L = 4 and Fig. 2b the L = 20 
case. It is clear that for the smaller lattice the system is probing both minima of the 
free energy {X^ = and X+ = 4.55 here), while for the larger lattice the system remains 
localized in the X = X^ well. As the lattice size is increased the throat separating the 
two minima becomes less and less dense until eventually the system becomes 'trapped' 
within one well. The time-scale for its eventual escape is much larger than any time-scale 
of interest in this problem, with 'Tqqq ^ OO clS 

Finally, in Fig. 3, we show the equilibrium values of the fraction fo{t), f^^, and 
of the volume-averaged order parameter, (X)'''', for a = 0.065 and A = 0.020, as a 
function of time step, 6t. Note that using too large a time step compromises the stability 
of the simulations, tending to drive the system's equilibrium towards the phase-mixed 
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symmetric state {f^ — > 0.5, (X)^ — X^g^). The results detailed below are therefore 
obtained using a time step 5t — 0.1. 

B. Numerical Results 

Based on the above discussion, we choose L — AS, Sx — 1, St — 0.1, and a — 0.065 
in all simulations. The experiment then consists in measuring the fraction of the volume 
in the 0-phase as a function of time for several values of A. As this involves a stochastic 
noise, we must perform an ensemble average over many runs in order to obtain physically 
reasonable results. 

In Fig. 4 we show the evolution of the ensemble- averaged fraction fo{t) for several 
values of A. It is clear that for small enough values of A the system remains well-locahzed 
in the 0-phase with /^^ ~ 1, while for larger values the two phases become completely 
mixed, with f^^ — > 0.5. Remarkably, the transition region between the two regimes is 
quite narrow, centered around A ~ 0.025. This can be seen from Fig. 5 where we show 
fo{t) for A = 0.024, 0.025, and 0.026. [The curves are noisier due to the fact that we 
must run for longer times in order to approach the equilibrium values of fo{t), being thus 
constrained to perform an ensemble average with fewer runs.] Note that for A = 0.026, 
f^^ ~ 0.5, while for A = 0.024, f^^ ~ 0.72. There is a pronounced change in the behavior 
of the system for A ~ 0.025. Furthermore, we find that the numerical curves can be fitted 
at all times by a stretched exponential, 

Ut) = (1 - f^) exp [-(VTeq)1 + f^"" , (22) 

where f^"^ is the final equilibrium fraction and Teq is the equilibration time-scale. In Table 
1 we hst a and Teq for several values of A. Note that for A = 0.025 the fit is obtained at 
late times by a power law, (smooth dashed curve in Fig. 5) 
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fo{t) |a=o.o25 oc t-' , (23) 

with k = 0.10(±0.02). This slowing down of the approach to equihbrium is typical of 
systems in the neighborhood of a second order phase transition, being known as 'critical 



slowing down' p4|. This behavior is suggestive of a phase transition between two possible 
regimes for the system, one in which the system is well-localized in the 0-phase, and the 
other in which there is a complete mixing between the two phases. Let us call these two 
regimes the 'strong' and the 'weak' regimes, respectively. Before we explore this idea 
any further, it is prudent to check if the final equilibrium fractions are sensitive to the 
viscosity parameter 77, which reflects the coupling of the system to the thermal bath. In 
Fig. 6 we show the approach to equilibrium for several values of 77 and for A = 0.020, 
which lies within the 'strong' regime. It is clear that whereas the equilibration time- 
scales are sensitive to the value of rj, with larger rj implying slower equilibration, the final 
equilibrium fractions are the same. This is precisely what one expects, as the coupling 
to the bath should not influence the final equilibrium properties of the system, but only 
how fast it equilibrates. 



Armed with these results, and invoking also the results in (2 + l)-dimensions [jT^, we 
define the equilibrium fractional population difference 

AF^^ie,) = ro^-fl'^ . (24) 

In Fig. 7 we show the behavior of AF'^^ as a function of A. There is a clear qualitative 
analogy between the behavior of AF'^^ as a function of A and the behavior of the magne- 
tization as a function of temperature in Ising models. Here, the order parameter is the 
equilibrium fractional population difference and the control parameter is A. Ac — 0.025 
is the critical value for the parameter A, which determines the degree of mixing of the 
system at Tc. 

We stress that the idea here is to probe the assumption of localization within the 
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0- phase as the system cools to Tc- Our results show that if the time-scales for cooling 
are slower than the equilibration time-scales of the system, for A > Ac there will be 
considerable phase mixing before the temperature drops below T^. This result can be 
made quite transparent by comparing the equilibrium value of the volume-averaged field 
(X)^'^, and the location of both the inflection point and the maximum of the potential 
with varying A. As can be seen from Fig. 8, the narrow transition region is clearly 
delimited by 

^inf < {Xy" <^max , (25) 

where X^-ai and are the inflection point and the maximum of the potential barrier, 
respectively. Note that for A > 0.026, f^"^ = 0.5 and (X)^i = X^^. RecaUing the 
information from Fig. 7, we conclude that there is a clear distinction between the 'strong' 
and 'weak' regimes. The existence of a potential barrier between the two phases at Tc 
does not imply in the system being localized in the 0-phase, for large enough values 
of A. This efficient thermal mixing of the phases will affect the dynamics of the phase 
transition as the temperature drops below Tc. For A > 0.025, if the cooling is slow 
enough, there is no reason to expect that the system will approach its final equilibrium by 
nucleation of bubbles larger than critical. Rather, the mechanism will resemble spinodal 
decomposition, where domains of the two phases will compete for dominance, with the 
-I— phase eventually dominating the volume. We stress that the value Ac = 0.025 is a weak 
lower bound; even if the system starts localized in the 0-phase at Tc for a smaller value of 
A, the potential barrier will also decrease for lower temperatures, and we should expect 
departures from nucleation settling in even for A < 0.025. [Of course, the decrease in 
temperature will also suppress thermal fluctuations.] The evaluation of the exact value of 
A for which nucleation theory will break down depends on the asymmetry of the potential, 
the temperature, and on the cooling rate, requiring further study. However, we can assert 
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that this value will be at least lower than Ac- 



V. CONCLUSIONS AND OUTLOOK 

In this paper we investigated the possibility that thermal fluctuations may induce con- 
siderable phase mixing for systems which exhibit two degenerate phases at their critical 
temperature. By modelling the nonequilibrium dynamics by means of a Langevin equa- 
tion with the system initially localized in one phase, we showed that complete mixing of 
the two phases can occur, despite the presence of a potential barrier between the phases. 
These results are of importance in the context of cosmological phase transitions, in par- 
ticular when the cooling rate is sufficiently slow compared to the equilibration time-scales 
of the system. This limit is implicit in our simulations, as we held the temperature fixed 
at Tc while the system equilibrated. Our results should also be of importance for systems 
studied in the laboratory. As we mentioned in the Introduction, binary fiuid mixtures 
separate by spinodal decomposition if the system is initially quenched with concentrations 
above the spinodal. The dynamics in that case is much slower though, and the transition 
between nucleation and spinodal decomposition is smoother [ITsf. However, we believe 
the essential physics to be the same, as we found that considerable phase mixing occured 
precisely as the equilibrium value of the volume-averaged scalar field (the equivalent of 
the concentration in the binary fiuid case) traversed within the 'transitional' region de- 
limited by the locations of the infiection point and the maximum of the potential barrier. 
Also, as mentioned in previous publications [1^, this phase mixing is characteristic of 



pre-transitional phenomena found in the isotropic-nematic transition of liquid crystals. 



which is known to be 'weakly' first order [25 



Inspired by the effective potential of the electroweak model, we measured the degree of 
mixing with respect to the value of the quartic self-coupling of the scalar order parameter. 
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However, as mentioned before, this is not a simulation of the electroweak transition, as 
our order parameter is a real scalar field. In fact, the critical value above which we found 
that the two phases mix, Ac ~ 0.025, is below the value of A which corresponds to a 
physical Higgs mass. Using the equivalence relation between the two models obtained in 
equation 16, a Higgs mass of Mh = 50 GeV would correspond to A = 0.0518, which is 
well within the 'weak' regime. For this value of A, the equilibrium value of the volume 
averaged field (X)'^^ = Xmax, and the two phases are indistinguishable. 

The results presented here are in qualitative agreement with a previous work in which 



phase mixing was investigated in (2+l)-dimensions ||T^. As in that case, we suggest that 
the distinction between 'weak' and 'strong' first order transitions be quantified in terms of 
the equilibrium fractional population difference between the two phases, AF'^'^ = fo^-~f+^- 
Within the limits of a finite lattice, we observed a sharp change in the behavior of AF"^"^ 
with respect to A, (see Fig. 7) which we suggest characterizes a second order phase 
transition between the 'symmetric' ('weak') phase AF^^ = and the 'broken-symmetric' 
('strong') phase AF'^'^ = 1. In order to study the details of the transition, a more thorough 
analysis of finite-size effects in the determination of the critical value Ac — 0.025 and of 
the critical exponent f3, obtained from the relation AF"^"^ = (Ac — A)^ in the neighborhood 
of Ac, must be performed. But taken together, the sharp change in AF^^ and the presence 
of critical slowing down about Ac provide substantial evidence for the existence of this 
transition. 

Finally, our work also calls for a more detailed analysis of the role of noise with 
spatio-temporal memory in simulations of the dynamics of phase transitions in field 
theories. The Markovian Langevin equation used here is certainly an approximation to 
more complicated couplings between the system and the thermal bath. It remains to be 
seen what role will colored noise play in the nonequilibrium dynamics of field theories. It is 
however hard to imagine that the nature of the noise can affect final equilibrium properties 
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of the system, although it may affect the equihbration time-scales. In this connection, 
we note that recent results by one of us show that at high enough temperatures, colored 
noise becomes white ||T9[ . Given our present level of understanding of the nonequilibrium 
dynamics of field theories it is only fair to expect many surprises in the coming years. 
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List of Figures 

Figure 1. The equilibrium values of the 0-phase fraction and of the volume- averaged 
field (X)^ are shown as a function of lattice length L for A = 0.020. 

Figure 2a. Phase space portrait for an arbitrarily chosen point of a lattice of length 
L — A. X — 0.010 in this run. Note how the system probes both minima of the potential. 

Figure 2b. Phase space portrait for an arbitrarily chosen point of a lattice of length 
L = 20. A = 0.010 in this run. Note how the system remains constrained within the 
0-phase. 

Figure 3. The equilibrium values of the 0-phase fraction and of the volume- averaged 
field (X)'^'^ are shown as a function of time step 6t for A = 0.020. 

Figure 4. Evolution of the ensemble-averaged fraction fo{t) for several values of A. 

Figure 5. Evolution of the ensemble-averaged fraction fo{t) for a narrow range about 
the critical value A = 0.025. The smooth dashed curve is the power law fit described in 
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equation 



Figure 6. Evolution of the ensemble-averaged fraction fo{t) for several values of the 
viscosity parameter 77 and for A = 0.020. Note how the final equilibrium fraction is inde- 
pendent of the value of 7]. 

Figure 7. The fractional equilibrium population difference AF^^ for several values of A. 
Note the sharp discontinuity about the critical value Ac — 0.025. 

Figure 8. Comparison of the equilibrium value of the volume-averaged field (X)*"' (cir- 
cles), with the location of the nearest inflection point to X = 0, Xi^f, and the poten- 
tial barrier, Xmaxj as a function of A. Note the existence of a transitional region for 
0.021 ~ A < 0.026. For A > 0.026, (X)'=^ = X^nax- 
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Table 1. The values of the equilibration time-scales and the exponents for the exponential 



fit of equation ^ for several values of A. 
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TABLE 1 



A 




a 


0.015 


25.0 


1.0 


0.020 


45.0 


1.0 


0.022 


60.0 


1.0 


0.024 


110.0 


0.80 


0.026 


220.0 


0.50 


0.028 


120.0 


0.75 


0.030 


100.0 


0.80 


0.035 


55.0 


0.80 


0.010 


10.0 


0.90 
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